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^ ! Abstract. We report extensive numerical simulations in the glass region 

for a simple model of short-ranged attractive colloids, the square well model. 
We investigate the behavior of the density autocorrelation function and of 
the static structure factor in the region of temperatures and packing fractions 
where a glass-glass transition is expected according to theoretical predictions. 
O ■ We strengthen our observations by studying both waiting time and history 

dependence of the numerical results. We provide evidence supporting the 
■r::j- I possibility that activated bond-breaking processes destabilize the attractive 

^ ■ glass, preventing the full observation of a sharp glass-glass kinetic transition. 



O ■ 1. Introduction 

O ■ 

The combination of theory [1, 2, 3, 4, 5, 6, 7], experiments [8, 9, 10, 11, 12] and nu- 
/\ ' merical simulations [13, 14, 15, 16, 18] for short-ranged attractive colloidal systems, i.e. 

c^ ■ colloid-polymer mixtures in which the size of the polymers is only a few percent that of 

the colloidal particles, has led in the last few years to a clear evidence of anomalous dy- 
namical properties for dense colloidal systems. In particular, the existence of a reentrant 
region of (super-cooled) liquid state, which is stable up to packing fractions higher than 
the correspondent purely repulsive system (hard spheres), has been established. From 
this liquid region one can approach dynamical arrest in two distinct ways [19], depend- 
ing on which control parameters for the system are manipulated. Thus, an increase of 
colloid density leads to the well-known hard-sphere (or simply 'repulsive') glass, while a 
stronger attraction, which is experimentally realized via an increase of polymer concen- 
tration, produces the so-called attractive glass, where colloidal particles are stabilized by 
'bonding'. Moreover, an efficient competition between these two mechanisms of arrest 
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originates a peculiar behaviour of dynamical correlation functions within the reentrant 
fluid region, in particular a logarithmic decay for the density autocorrelation functions 
and a sub-diffusive regime for the mean squared particle displacement [15, 16, 27]. All 
these features have been firstly predicted by the Mode Coupling Theory (MCT) [21] on 
the basis of the existence in the control parameter space for these systems of a higher or- 
der singularity in the solutions of MCT equations, which entirely governs the anomalous 
dynamics. This particular point in the phase diagram, named A4 singularity, is found for 
a particular value of the three control parameters, i.e. temperature, volume fraction of 
the colloids and range of the attractive interaction, where the liquid, the attractive glass 
and the hard-sphere glass solutions all merge into a single point. The key parameter for 
the existence of such a point is the attractive range of the potential A (experimentally 
set by the ratio of polymer gyration radius to the colloidal radius): for A > A* there is 
no distinction between the two glassy states, while for A < A* there exists a discontin- 
uous glass-glass transition between the two [3, 4]. The glass-glass transition terminates 
at an end-point, named A^ beyond which the two glasses are indistinguishable. Such 
kinetic glass-glass transition is one of the most difficult MCT prediction to detect, due 
to the fact that it lies entirely in a non-equilibrium region where both experiments and 
simulations are highly non-trivial. 

On the other hand, while there is convincing evidence, also for molecular glass- 
formers, of the accuracy of MCT predictions in the liquid regime[22], the same does not 
apply for regions where the theory predicts a glass. Indeed, the MCT glass is commonly 
termed 'ideal', because it arises without the inclusion of the so-called hopping processes, 
i.e. further relaxing processes which eventually allow the system to still restore ergodicity, 
beyond the ideal MCT transition [23]. This may be better explained in the free energy 
landscape picture, where a glassy state is represented as a local (metastable) minima 
and a hopping process corresponds to a jump between different minima, via a low order 
saddle-point [24]. 

In this paper, we report an extensive numerical study of the glass-glass transition 
for a specific simple model [15], for which the location of the MCT transition lines and of 
the higher order singularity have been extensively studied and estimated via a mapping 
procedure [16, 27]. Thereby, after explaining in detail how this mapping is performed, 
we report results in the glassy region, including static and dynamical properties of the 
glassy states, as well as we compare results obtained following different histories for 
the system, due to their non-equilibrium nature. This will conffim the importance of 
hopping processes in some mechanisms of glassification, and provide more insight on the 
MCT applicability to describe different mechanism of arrest. A short account of the 
main results was already anticipated in a recent Letter [29]. 

2. Methods 

We perform event-driven Molecular Dynamics simulations for a binary mixture of parti- 
cles interacting via a narrow square-well (SW) potential. The parameters of the mixture 
are chosen in order to avoid crystallization in the system. Thus, we consider a 50:50 
mixture of 700 particles of mass m with diameters aAA = 1-2 and asB = 1; the hard- 
core diameter for the AB interaction is cxab = (o'aa + o'bb)/2. We fix the width of 
the attractive well as Aij/[a)ij + Aj^) = 3%. Temperature is measured in units of the 
well-depth uq, while time units are (Tbb(^/mo)^^^- 



This system was extensively studied previously by means of theory (MCT) and 
simulations. The combination of both tools allowed us to determine, within some ap- 
proximations, the location of the MCT lines, and consequently of MCT higher order 
singularities, for this particular system. In this way, we obtain the best possible diagram 
for studying the dynamics of the system near the A4 singularity [16] and crossing the 
glass-glass transition [29]. 

Here, we briefly review how this mapping procedure between theory and simula- 
tions is performed. Firstly, we solved MCT equations for our particular binary mixture 
using as inputs the partial static structure factors calculated within Percus-Yevick (PY) 
approximation [28] . In particular, for the static structure calculations, we solved numeri- 
cally the Ornstein-Zernike equation on a grid of 20000 wavevectors, with mesh 0.3141593, 
while for the binary MCT equations we used a grid of 2000 wavevectors with the same 
mesh. By solving the equations for a large range of T — values, the PY-MCT glass lines 
have been detemined. Next, we calculated the diffusion coefficients from simulations in 
the liquid region, covering 4 decades in diffusivity values. A fit of these values according 
to the power-law D ~ (0 — (j)g[T)y, provides estimates for the glass transition packing 
fraction (/)g{T) (at which the diffusion coefficient would vanish) for the various isotherms 
and a characteristic exponent '^ which, according to MCT predictions, should grow near 
a higher order singularity. Note that 7 is related to the exponent parameter A which 
becomes exactly equal to 1 at the MCT singularities (^3,^4). Figure 1 shows the fit of 
the diffusivity data, together with the parameters 7 and A obtained from them. Errors 
in these values are given by the errors in the fits [30]. The locus 0g(T) provides the best 
numerical estimate of the glass transition lines. Next we perform a bilinear transfor- 
mation in the control parameter space (0,T), to superimpose the PY-MCT results for 
the ideal glass line with the extrapolated results from the simulations, following an idea 
proposed by Sperl [27]. The resulting transformation is found to be, 

0,,™^ 1.897 0MCT- 0.3922; 

Ts,m -^ 0.5882 Tmct - 0.225. (1) 

This transformation provides a tool for comparing theoretical predictions and simula- 
tion results. Improvements over such procedure could be obtained by using the numerical 
"exact" structure factors as oppose to the PY-closure in the estimate of the ideal MCT 
lines. Indeed, despite the quite good agreement of the PY solutions, small differences 
between PY and numerical structure factors are crucial for a full quantitative agreement 
of theory and simulations results, as recently discussed in [31]. Assuming that the bilin- 
ear transformation is invariant for small changes of A, we can locate theoretically the 
position of the A4 singularity, moving along the A-axis in the control parameter space 
and find the corresponding location for the real system [16]. The qualitative general 
agreement between predictions and simulations results for the mean squared displace- 
ment and density auto-correlation functions shown in Ref. [16] strongly supports the 
validity of the mapping procedure. Having established the quality of the mapping, we 
can locate the theoretically predicted glass-glass transition for the A = 0.03 case which 
refer to the present simulation parameters. More specifically, the glass-glass line is found 
to be located in the (0, T)-plane approximately between (0.625,0.37) and (0.64,0.41). 
The resulting liquid-glass and glass-glass lines for the model are shown in Fig. 2. 

Our aim is to perform MD simulations across the glass-glass transition line and 
study the dynamical properties, in particular the collective density-density correlation 



functions, both in the hard-sphere and in the attractive glass region. A discontinuity in 
such properties should be observed, according to MCT, at the glass-glass transition. In 
order to access the glass region, we used the following procedure, that is schematized in 
Fig. 2. We consider the slowest equilibrated liquid state point in the reentrant region, i.e. 
corresponding to initial packing fraction and temperature of (pi = 0.612 and Ti/uo = 0.6 
(labelled as A in the figure). We then rapidly compress the configuration by progressively 
growing the particle radii in successive steps. In order to investigate the glass-glass 
transition we choose to compress isothermally the system up to = 0.625 and up to 
4> = 0.635. At each of the two packing fractions, we then suddenly change the mean 
square velocities to achieve configurations at several different temperatures (indicated 
with crosses). A thermostat with a very short characteristic time kept the temperature 
constant during the following aging evolution. This procedure has been repeated for 136 
independent initial configurations extracted from the the initial equilibrium liquid state. 
Being the system out-of-equilibrium, the dynamical observables will have a depen- 
dence also on the waiting time t^, and on the chosen initial state [25, 26]. Thus, in 
particular the density autocorrelation functions are defined as, 

0,(t»+t,t^) = (P;(t^+t)pq(t^))/(|pq(t»)|') (2) 

where Pq(t) = -j^J2i^^'^^'^^^ and rj(t) are the coordinates of particle i at time t. The 
correlation functions calculated in this paper always refer to all particles in the mixture. 
Averages are taken over up to 136 independent initial configurations from the initial 
liquid state. To assess the problem of a possible path-dependence of the results, we 
also consider a different initial liquid state, this time within the attractive region, i.e. 
(pi = 0.58, Ti/uQ = 0.3 (labelled as B in the figure). Moreover, we check that our results 
are waiting-time independent within a well-defined, even though narrow, time window. 
This was already discussed in [29]. For the present purposes, we recall that the results 
reported in this article correspond to t^ = 4 ■ 10^ which provides sufficient equilibration 
for the decay processes taking place on a time scale one order of magnitude smaller [26]. 
For this reason we only show the decay of correlation up to t = 400. However, in order to 
provide convincing evidence that results for t < 400 do not depend on the chosen t^, we 
will also report results corresponding to a t^j larger by an order of magnitude. Finally, 
we will discuss the behaviour of static structure factors in the glassy region. 



3. Results 

Before showing the simulation results, we start by reporting MCT predictions. In Figure 
3, we show the behaviour of the (total) non-ergodicity factor fg = limt^oo (pqit) (where we 
assume no waiting time dependence), crossing along an isochore the glass-glass transition, 
at temperatures whose distance from the transition is the same as for the simulated 
points. It is evident that fq does not change very much in the hard-sphere glass region, 
but jumps abruptly at the glass-glass transition, and then changes dramatically with 
temperature in the attractive region. We here want to draw the reader's attention to 
the non-monotonic behaviour of fq in the hard-sphere region with temperature, as this 
non-trivial feature will be found also in the simulations. The crossing arises because at 
this isochore at intermediate temperatures (see Fig. 2), there are two distinct transitions 
close-by, i.e. standard liquid-glass and glass-glass one, which can effectively compete. 



Similarly, if we plot the full time evolution of (f)q{t) at the same points for a particular 
wave-vector, eg. at the static structure factor first peak, we find a net difference in the 
levels of the reached plateaux, i.e. very high for the attractive glasses and much lower 
for the hard-sphere ones (see Fig. 2b of ref. [29]). We recall here that the glass-glass 
transition is a purely kinetic phenomenon, which arises from the non-linearity of MCT 
equations. Thus, there is no abrupt change in the static structure, passing from one 
glass to the other. 

Now we turn to the simulation results. We plot the density autocorrelation func- 
tions (pq(t, t^) for = 0.625 at all studied temperatures for the representative wave- vector 
qasB — 25, since their behaviour is similar for all studied wave-vectors. The following 
features are evident in Figure 4: 

• at higher temperatures, the system is already able to relax to its plateu, even within 
the small studied time window. Also, the predicted non-monotonic behaviour of 
the plateau with temperature, found theoretically, is reproduced in full details; 

• lowering the temperature, there is a smooth change in the correlators, which start 
to display a quite evident logarithmic decay, signature of the higher order MCT 
singularity that is closeby; 

• at temperatures already below the glass-glass line, i.e. for T/uo < 0.4, the cor- 
relators continue quickly to drop down, below the expected attractive plateau; a 
further decrease in temperature elongates the time duration of the plateau, but a 
decay at longer times is present even in the extreme case T/uq = 0.1; 

• looking in more details, we observe that theoretically fq{T /uq = 0.9) < fq{T /uq = 
0.65) (see Fig. 3) while the simulations show the opposite feature at our final time of 
observation. This could be probably due to the inaccuracy of PY structure factors 
that we used to solve MCT equations [31], or to the fact that hopping processes 
start already to produce deviations from MCT at T/uq = 0.65 (in this respect, at 
the immediately lower temperature T/uq = 0.5, MCT predictions are completely 
far off the simulation data). 

A very similar behaviour is found analyzing the different packing fraction = 0.635, 
already reported in [29]. The only difference with the previous case is that no evident 
crossing of fg was found there. This is again in agreement with MCT, because the 
influence of the liquid-glass transition is now less strong, since the chosen isochore lies 
more deeply into the glassy region. 

Conclusions that can be drawn from these results are that the hard-sphere glass 
predictions are fully reproduced by simulations, while the attractive glass seems to be 
less stable than predicted. An interpretation of these findings is possible taking into 
account a relevant parameter for the system at low temperatures, i.e. the lifetime of the 
attractive bonds. The square well model is particularly suited to study bond dynamics, 
since a bond is unambigously defined when two particles have a pair interaction energy 
equal to —Uq, i.e. they are at a distance less or equal than A. To monitor the dynamical 
evolution of the bonds, we introduce the 'bond' correlation function (psitw+t, t^), defined 
as 

(pB(ty, + t,ty,) = (^nij(tw)nij(ty, + t)) / [N u B (t u,)] (3) 

i<j 



where nij{tw) is 1 if two particles are bonded and otherwise, while nB{tw) = 
{I2i<j''^ij(tw))/N is the average number of bonds per particle at tw (pB counts how 
many of the bonds found at time t^ are still present after time t^ + t, independently 
from any breaking-reforming intermediate process. 

The time evolution of (psitw + t^tw), shown in Fig. 5, is very similar to the one of 
4'q{'^w + t, tw), suggesting that the decay of density- density correlation is associated to the 
dynamics of bond breaking. The reason why the shape of the correlators in Fig. 5 looks 
similar to density correlators even for high temperatures, lies in the fact that we have 
taken into account also bonds re-forming, once they are broken. Indeed, the packing 
fraction we are studying (0 = 0.635) is well above the glass transition packing fraction 
{(f)g ~ 0.58) for simple hard spheres, thus particles are confined to rattle in their cages, 
and bonds gets continuously broken and reformed. To give a convincing evidence of the 
different dynamics taking place at high and low temperatures, we report in Fig. 6 the 
bond correlation function, at the two extreme temperatures T/uq = 0.1 and T/uq = 1.5, 
defined above, as well as the same functions but where nij{tw) is set to zero at the time 
when the ij-bond is broken, and kept to zero from then on, irrespectively of the presence 
of the ij-bond at later times. In the high temperature case, the bond correlation function 
goes almost to zero when broken bonds are not re-counted, while at low temperature 
there is no substantial change. In our previous work [29], we have already put forward 
the observation of relevance of bond dynamics for attractive glasses, by showing how a 
characteristic type of hopping processes takes place in the attractive glass, preempting 
the observation of a sharp glass-glass transition, as predicted by MCT, and destabilizing 
the attractive glass. The study of the effect of the bonds lifetime on the full dynamical 
behaviour of short-ranged attractive colloids has also been explored in [32] . 

In the present work, we support the hypothesis of bond-breaking hopping processes 
preempting the observation of a sharp glass-glass transition, as predicted by MCT, show- 
ing that results for the density correlators as shown in Fig. 4 do not depend on the chosen 
waiting time or on the history path dependence. 

Thereby, we first report in Fig. 8 the density correlation functions for three repre- 
sentative temperatures, i.e. T = 1.5, 0.35, 0.1, for a waiting time of 5 ■ 10^, compared to 
the case examined so far for a tu, smaller by a full order of magnitude. These results are 
averaged over up to 40 configurations only, due to the very large computer time involved 
in the calculation (approximately 10^ CPU hours for each temperature). It is evident 
that no significant change in the density correlators behaviour arises, though small dis- 
crepancies are present also due to the larger errors present in the longer simulation runs. 
However, the most important feature to underline is that an increase in waiting time 
does not manifest in an increase of the time duration of the attractive plateau. Our 
study shows that the stability of such plateau, i.e. of the bond lifetime, can only be 
achieved by decreasing temperature. Thus, the observed phenomenon is entirely due to 
bond dynamics, which is independent on the waiting time and completely controls the 
attractive glass transition. 

Of course, an interesting issue would be to know exactly where, at infinite times, 
the attractive glass correlation functions would decay to. In other words, we know that 
the plateau must be finite, and its lower bound should be represented by the hard-sphere 
plateau. So, we can speculate that the long-time cage of the attractive glass is still of 
the order of the hard-sphere one, though at each time particles are trapped within the 
attractive well-width, which is highly affected by breaking processes, differently from 
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simple packing. Unfortunately, we cannot support this picture by simulations, since it 
would require many additional decades in time, beyond computational efforts. 

The last point we wish to address is the dependence on the chosen liquid ini- 
tial state. Indeed, one might argue that we started our non-equilibrium (compres- 
sion-|-quenching) path from a state point close to a hard-sphere glass. Thus, we can 
repeat the full procedure from a different initial state close to the attractive glass tran- 
sition. For this purpose, we refer to the initial state point B of Figure 2. Results for 
the density correlators are compared, for the two different initial points, in Figure 8, for 
temperatures T/uo = 0.1, 0.4, 0.9. It is clear that quantitatively results are different, but 
qualitatively they are not. As expected, starting from a lower temperature causes the 
produced glass to have a higher correlation function, i.e. it imprints in it more of the 
attractive features. However, after the microscopic time, the decay is very similar in the 
two cases. Most importantly, the lowest temperature plateau is not stabilized by the 
different path. 

It is interesting to note that if one examines the complete T-dependence of the 
density correlators at the same isochore = 0.625 as in Fig. 4, but now starting from B, 
one does not observe the same non-monotonic behaviour, at high temperatures, for the 
density correlators. This is due to the fact that we are starting from a configuration close 
to the attractive line, which somehow remains imprinted in the memory of the system 
(see Figs. 8,9). Thus, probably, our small time window of observation is not sufficient 
for the system to forget about its initial state, which prevents the observation of smaller 
plateaux especially at higher temperatures. In other words, the system cannot really feel 
enough competition from two transitions, and the only relevant one for the dynamics 
here is the glass-glass one. Interestingly enough, this behaviour resembles closely the 
one at the highest isochore (0 = 0.635) starting from A[29]. 

Finally, we want to focus on the behaviour of the static structure factors in the 
glassy region. As explained above, S{q) is the input needed to solve MCT equations, 
and at the glass-glass transition it does not have any peculiar change, other than the 
smooth one characteristic of a slight change of thermodynamic control parameters. In 
Figure 9, we show the (total) static structure factors for T/uq = 0.1,0.9 together with 
the equilibrium configuration one, both for the path starting from point A (top panel) 
and from point B (bottom panel). It is clear that the structure just changes smoothly 
along the quenching paths, but does change with respect to the initial equilibrium state 
because of the density increase. In the glass, no detectable waiting time dependence for 
the structure is observed. Static structure factors between the two panels, corresponding 
to the same temperatures after different paths, cannot be really distinguished within 
numerical errors. However, we expect small changes due to the different imprinted 
structures for the initial state-point. This could help explaining why the results for the 
dynamical correlators are different between the two paths (see Fig. 8.). The results 
shown here for the static structure factors are quite different than those found in a 
micellar system [11], and that could explain the different conclusions drawn in that 
paper. 

4. Conclusions 

In this paper, we have reported numerical simulations along the theoretically predicted 
glass-glass transition between repulsive and attractive glass for a simple model of short- 



ranged attractive colloids. A simple mapping between theory and simulations, performed 
for this particular model, allowed us to have an accurate location of this glass-glass 
transition, and thus, to have a consistent test of the theoretical predictions by simulations 
in the glass. 

Studying the glassy dynamics, we can discriminate time intervals where waiting 
time does not play a significant role, and where there is only a minor dependence on the 
preparation of the glassy state. It is, therefore, an ideal system to study and compare 
the two mechanisms of glassification, i.e. packing and attraction. 

Our results clearly show that a glass-glass transition exists. Indeed, there is a 
significant difference in the dynamical behaviour between hard-sphere and attractive 
glass, though the crossover from one kind to the other is smoother than predicted by 
MCT. Thus, the transition does not manifest as sharply as predicted by the theory, due 
to bond-breaking processes that reduce the stability of the attractive glass. The detailed 
studies of waiting time and path dependence strengthen this interpretation, because they 
do not have any effect on the stability of the attractive plateau. This, indeed, can only 
be made longer- lasting in time by a decrease of temperature, i.e. an (Arrhenius) increase 
of the bond lifetime. Therefore, we claim that bond-breaking activated processes have 
the same relevance on attractive glass as they do in molecular glass-formers, while hard- 
sphere glasses, due to their entropic origin, are basically unaffected by them. On the other 
hand, the intrinsic different nature of the two glasses can still be experimentally detected. 
In particular, for systems where the bond lifetime is long enough, or equivalently by 
looking at time-scales that are shorter than the bond lifetime, the MCT predictions for 
the glass-glass transition could be fully recovered. 
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Figure 1. Left panel: Diffusion coefficients in the liquid region for variuos 
studied temperatures. Dashed lines are power-law fits, with parameters (pg, 
i.e. the glass transition packing fraction at which D = 0, and the expo- 
nent 7. Right panel: behaviour of 7, and of the exponent parameter A with 
temperature. 
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Figure 2. Map of the simulations. We consider two initial state points, la- 
beled A and B, lying along the iso-diffusivity curve (diamond - small dashed 
line) calculated in [15], respectively close to the hard-sphere and to the attrac- 
tive glass transition. We compress these configurations up to two isochores, 
= 0.625 and = 0.635, and follow quenches at various temperatures (crosses 
- dotted lines). The solid curve is the mapped glass line from Ref. [16], as 
described in the text. 
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Figure 3. MCT predictions for tlie number density non-ergodicity parameter 
fq for ~ 0.5362 (corresponding to (psiM = 0.625). The top one refers to 
high temperatures, in the hard-sphere glass region, and the bottom one to 
the lower ones in the attractive glass. The case T/uo = 0.4 is reported in 
both panels. Note the non-monotonicity of fg in the hard-sphere region, as 
well as the difference of one order of magnitude in the qa scale betwen the 
two panels. 
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Figure 4. Density correlators calculated from the simulation at 

for qaBB ^ 25 and ^ = 3754. 
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Figure 5. Time evolution of the bond correlators (psitw + t, t^) for = 0.635 
and t^. Temperatures are the same as in the previous figure (corresponding 
to same symbols). 




Figure 6. Time evolution of the bond correlators 0B(t«, + t, tw) for = 0.635 
and t^j for temperatures T/uq = 1.5 and T/uq = 0.1, comparing the case 
when bonds after being broken are still counted (same as 5, dashed lines) or 
not (curves labelled as 'broken', full lines). 
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Figure 7. Study of the waiting time dependence for three different tempera- 
tures at = 0.635. Points refer to t^ = 50344, one order of magnitude larger 
than t^. Clearly, the plateaux are not affected by the difference in waiting 
time, at all temperatures. 
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Figure 8. Study of the history dependence for three different temperatures 
at = 0.635 and t^. Points refer to initial state point B in Fig. 2, while solid 
lines are for the standard initial state point A. 
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Figure 9. Static structure factors at different temperatures for = 0.635, 
compared to each other and to the equihbrium static structure factor of the 
initial hquid: top panel is for initial state point A {Ti = 0.6, 0j = 0.612), 
bottom for initial state point B (Ti = 0.3, (pi = 0.58). 
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